1.1.1 Workshop Materials

All materials for this workshop is available in my standard GitHub repo:

https://github.com/kaybenleroll/dublin_r_workshops

book cover

The content of this workshop is partly based on the book “Applied Survival Analysis Using R” by Dirk F. Moore. The data from this book is available from CRAN via the package asaur and there is a GitHub repo for the code in the book also:

https://github.com/kolaczyk/sand

1.1.2 Example Datasets

This workshop will use a number of different time-to-event datasets, so we look at those first.

telco_churn_tbl %>% glimpse()

## Observations: 5,000
## Variables: 21
## $ state            <chr> "KS", "OH", "NJ", "OH", "OK", "AL", "MA", "MO", "L...
## $ accountdur       <dbl> 128, 107, 137, 84, 75, 118, 121, 147, 117, 141, 65...
## $ areacode         <chr> "415", "415", "415", "408", "415", "510", "510", "...
## $ phonenumber      <chr> "382-4657", "371-7191", "358-1921", "375-9999", "3...
## $ intlplan         <chr> "no", "no", "no", "yes", "yes", "yes", "no", "yes"...
## $ vmailplan        <chr> "yes", "yes", "no", "no", "no", "no", "yes", "no",...
## $ vmailmessage     <dbl> 25, 26, 0, 0, 0, 0, 24, 0, 0, 37, 0, 0, 0, 0, 0, 0...
## $ daymins          <dbl> 265.1, 161.6, 243.4, 299.4, 166.7, 223.4, 218.2, 1...
## $ daycalls         <dbl> 110, 123, 114, 71, 113, 98, 88, 79, 97, 84, 137, 1...
## $ daycharge        <dbl> 45.07, 27.47, 41.38, 50.90, 28.34, 37.98, 37.09, 2...
## $ evemins          <dbl> 197.4, 195.5, 121.2, 61.9, 148.3, 220.6, 348.5, 10...
## $ evecalls         <dbl> 99, 103, 110, 88, 122, 101, 108, 94, 80, 111, 83, ...
## $ evecharge        <dbl> 16.78, 16.62, 10.30, 5.26, 12.61, 18.75, 29.62, 8....
## $ nightmins        <dbl> 244.7, 254.4, 162.6, 196.9, 186.9, 203.9, 212.6, 2...
## $ nightcalls       <dbl> 91, 103, 104, 89, 121, 118, 118, 96, 90, 97, 111, ...
## $ nightcharge      <dbl> 11.01, 11.45, 7.32, 8.86, 8.41, 9.18, 9.57, 9.53, ...
## $ intlmins         <dbl> 10.0, 13.7, 12.2, 6.6, 10.1, 6.3, 7.5, 7.1, 8.7, 1...
## $ intlcalls        <dbl> 3, 3, 5, 7, 3, 6, 7, 6, 4, 5, 6, 5, 2, 5, 6, 9, 4,...
## $ intlcharge       <dbl> 2.70, 3.70, 3.29, 1.78, 2.73, 1.70, 2.03, 1.92, 2....
## $ custservicecalls <dbl> 1, 1, 0, 2, 3, 0, 3, 0, 1, 0, 4, 0, 1, 3, 4, 4, 1,...
## $ churned          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F...

telco_churn_tbl %>% head()

## # A tibble: 6 x 21
##   state accountdur areacode phonenumber intlplan vmailplan vmailmessage daymins
##   <chr>      <dbl> <chr>    <chr>       <chr>    <chr>            <dbl>   <dbl>
## 1 KS           128 415      382-4657    no       yes                 25    265.
## 2 OH           107 415      371-7191    no       yes                 26    162.
## 3 NJ           137 415      358-1921    no       no                   0    243.
## 4 OH            84 408      375-9999    yes      no                   0    299.
## 5 OK            75 415      330-6626    yes      no                   0    167.
## 6 AL           118 510      391-8027    yes      no                   0    223.
## # ... with 13 more variables: daycalls <dbl>, daycharge <dbl>, evemins <dbl>,
## #   evecalls <dbl>, evecharge <dbl>, nightmins <dbl>, nightcalls <dbl>,
## #   nightcharge <dbl>, intlmins <dbl>, intlcalls <dbl>, intlcharge <dbl>,
## #   custservicecalls <dbl>, churned <lgl>

prostate_surv_tbl %>% glimpse()

## Observations: 14,294
## Variables: 5
## $ grade    <fct> mode, mode, poor, mode, mode, poor, poor, mode, mode, poor...
## $ stage    <fct> T1c, T1ab, T1c, T2, T1c, T2, T1c, T1ab, T1c, T2, T1c, T1c,...
## $ ageGroup <fct> 80+, 75-79, 75-79, 70-74, 70-74, 75-79, 80+, 80+, 75-79, 7...
## $ survTime <int> 18, 23, 37, 27, 42, 38, 18, 78, 47, 13, 81, 23, 21, 13, 11...
## $ status   <int> 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...

prostate_surv_tbl %>% head()

## # A tibble: 6 x 5
##   grade stage ageGroup survTime status
##   <fct> <fct> <fct>       <int>  <int>
## 1 mode  T1c   80+            18      0
## 2 mode  T1ab  75-79          23      0
## 3 poor  T1c   75-79          37      0
## 4 mode  T2    70-74          27      0
## 5 mode  T1c   70-74          42      0
## 6 poor  T2    75-79          38      2

pharmaco_smoker_tbl %>% glimpse()

## Observations: 125
## Variables: 14
## $ id             <int> 21, 113, 39, 80, 87, 29, 16, 35, 54, 70, 84, 85, 25,...
## $ ttr            <int> 182, 14, 5, 16, 0, 182, 14, 77, 2, 0, 12, 182, 21, 3...
## $ relapse        <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0...
## $ grp            <fct> patchOnly, patchOnly, combination, combination, comb...
## $ age            <int> 36, 41, 25, 54, 45, 43, 66, 78, 40, 38, 64, 51, 37, ...
## $ gender         <fct> Male, Male, Female, Male, Male, Male, Male, Female, ...
## $ race           <fct> white, white, white, white, white, hispanic, black, ...
## $ employment     <fct> ft, other, other, ft, other, ft, pt, other, ft, ft, ...
## $ yearsSmoking   <int> 26, 27, 12, 39, 30, 30, 54, 56, 25, 23, 30, 35, 23, ...
## $ levelSmoking   <fct> heavy, heavy, heavy, heavy, heavy, heavy, heavy, lig...
## $ ageGroup2      <fct> 21-49, 21-49, 21-49, 50+, 21-49, 21-49, 50+, 50+, 21...
## $ ageGroup4      <fct> 35-49, 35-49, 21-34, 50-64, 35-49, 35-49, 65+, 65+, ...
## $ priorAttempts  <int> 0, 3, 3, 0, 0, 2, 0, 10, 4, 10, 12, 1, 5, 6, 5, 2, 1...
## $ longestNoSmoke <int> 0, 90, 21, 0, 0, 1825, 0, 15, 7, 90, 365, 7, 1095, 1...

pharmaco_smoker_tbl %>% head()

## # A tibble: 6 x 14
##      id   ttr relapse grp     age gender race  employment yearsSmoking
##   <int> <int>   <int> <fct> <int> <fct>  <fct> <fct>             <int>
## 1    21   182       0 patc…    36 Male   white ft                   26
## 2   113    14       1 patc…    41 Male   white other                27
## 3    39     5       1 comb…    25 Female white other                12
## 4    80    16       1 comb…    54 Male   white ft                   39
## 5    87     0       1 comb…    45 Male   white other                30
## 6    29   182       0 comb…    43 Male   hisp… ft                   30
## # ... with 5 more variables: levelSmoking <fct>, ageGroup2 <fct>,
## #   ageGroup4 <fct>, priorAttempts <int>, longestNoSmoke <int>

1.2.1 Data Censoring and Truncation

In the problems we work on, our data observation process does not allow us to fully observe the variable in question. Instead, our observations are often right-censored - that is, we know the value in question is at least the observed value, but it may also be higher.

Expressing this formally, suppose \(T^*\) is the time to event, and \(U\) is the time to the censoring, the our observed time \(T\), and censoring indicator, \(\delta\), are given by

\[\begin{eqnarray*} T &=& \min(T^*, U) \\ \delta &=& I[T^* \leq U] \end{eqnarray*}\]

A less common phenomenon is left-censoring - where we observe the event to be at most a given duration. This may happen in medical studies where continual observation of the patients is not possible or feasible.

1.2.2 Hazard and Survival Functions

Our key goal is to find the survival distribution - the distribution of times from some given start point to the time of the event.

Two common ways of specifying this distribution are the survival function, \(S(t)\) and the hazard function, \(h(t)\).

\(S(t)\) is the probability of surviving to time \(t\), so is defined as follows:

\[ S(t) = P(T > t), \, 0 \leq t \leq \infty \]

Thus, \(S(0) = 1\) and decreases monotonically with time \(t\). As it is a probability it is non-negative.

We can also define the survival function in terms of the instantaneous failure rate, the probability of failing at exactly time \(t\). More formally,

\[ \lambda(t) = \lim_{\delta \to 0} \frac{P(t \leq T \leq t + \delta \; | \, T > t)}{\delta} \]

This is also called the intensity function or the force of mortality

1.2.3 Cumulative Functions

Some analyses make it easier to use the cumulative hazard function

\[ \Lambda(t) = \int^t_0 \lambda(u) \, du \]

As a consequence of this definition we also have

\[ S(t) = \exp(- \Lambda(t)) \]

1.2.4 Mean and Median Survival

Finally, two common statistics for survival are the mean survival time and the median survival time:

The mean survival time is the expected value of the distribution, using standard probability definitions:

\[ \mu = E(T) = \int^\infty_0 t f(t) \, dt = \int^\infty_0 S(t) \, dt \]

The median survival time, \(t_{\text{med}}\) such that

\[ S(t_{\text{med}}) = 0.5 \]

1.3.1 Constant Hazard

t_seq <- 1:100

h_const_seq <- rep_along(t_seq, 0.01)
S_const_seq <- cumprod(c(1, (1 - h_const_seq)))


hazard_plot <- ggplot() +
    geom_line(aes(x = t_seq, y = h_const_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Instantaneous Hazard")

cuml_plot <- ggplot() +
    geom_line(aes(x = c(0, t_seq), y = S_const_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Cumulative Survival")

plot_grid(hazard_plot, cuml_plot, nrow = 2)

If we increase the constant hazard, we get a similar shape, but the curve declines faster.

h_const_high_seq <- rep_along(t_seq, 0.03)
S_const_high_seq <- cumprod(c(1, (1 - h_const_high_seq)))


hazard_plot <- ggplot() +
    geom_line(aes(x = t_seq, y = h_const_high_seq)) +
    geom_line(aes(x = t_seq, y = h_const_seq)
             ,colour = 'blue', linetype = 'dashed') +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Instantaneous Hazard")

cuml_plot <- ggplot() +
    geom_line(aes(x = c(0, t_seq), y = S_const_high_seq)) +
    geom_line(aes(x = c(0, t_seq), y = S_const_seq)
             ,colour = 'blue', linetype = 'dashed') +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Cumulative Survival")

plot_grid(hazard_plot, cuml_plot, nrow = 2)

1.3.2 Early Hazard

early_seq <- seq(0.10, 0.01, by = -0.01)
late_seq  <- rep(0.01, 100 - length(early_seq))

h_early_seq <- c(early_seq, late_seq)
S_early_seq <- cumprod(c(1, (1 - h_early_seq)))


hazard_plot <- ggplot() +
    geom_line(aes(x = t_seq, y = h_early_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Instantaneous Hazard")

cuml_plot <- ggplot() +
    geom_line(aes(x = c(0, t_seq), y = S_early_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Cumulative Survival")

plot_grid(hazard_plot, cuml_plot, nrow = 2)

1.3.3 Late Hazard

late_seq  <- seq(0.01, 0.10, by = 0.0015)
early_seq <- rep(0.01, 100 - length(late_seq))

h_late_seq <- c(early_seq, late_seq)
S_late_seq <- cumprod(c(1, (1 - h_late_seq)))


hazard_plot <- ggplot() +
    geom_line(aes(x = t_seq, y = h_late_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Instantaneous Hazard")

cuml_plot <- ggplot() +
    geom_line(aes(x = c(0, t_seq), y = S_late_seq)) +
    expand_limits(y = 0) +
    xlab("Time") +
    ylab("Cumulative Survival")

plot_grid(hazard_plot, cuml_plot, nrow = 2)

2.1.1 Using Kaplan-Meier

In R, we construct KM estimators using the survfit() function.

Before we move on to our datasets, we start with a small set of data.

tt   <- c(7, 6, 6, 5, 2, 4)
cens <- c(0, 1, 0, 0, 1, 1)
grp  <- c(0, 0, 0, 1, 1, 1)

Surv(tt, cens)

## [1] 7+ 6  6+ 5+ 2  4

sample_tbl <- data_frame(tt = tt, cens = cens, grp = grp)

example_km <- survfit(Surv(tt, cens) ~ 1, data = sample_tbl, conf.type = 'log-log')

plot(example_km)

Basic plotting routines are worth trying, but the survminer package has specialised plots that use ggplot2 to create them.

ggsurvplot(example_km)

Printing out the ‘fitted’ object gives us some basic statistics:

example_km %>% print()

## Call: survfit(formula = Surv(tt, cens) ~ 1, data = sample_tbl, conf.type = "log-log")
## 
##       n  events  median 0.95LCL 0.95UCL 
##       6       3       6       2      NA

We get more details from the summary() function:

example_km %>% summary()

## Call: survfit(formula = Surv(tt, cens) ~ 1, data = sample_tbl, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     2      6       1    0.833   0.152       0.2731        0.975
##     4      5       1    0.667   0.192       0.1946        0.904
##     6      3       1    0.444   0.222       0.0662        0.785

2.1.2 Follow-up Time

A useful measure may be how long the observation period lasts, something that can be subtly difficult to measure.

One method is to switch the censoring labels - that is, we consider the original event as a censoring of the observation in the study.

sample_tbl <- sample_tbl %>%
    mutate(follow = 1 - cens)

follow_km <- survfit(Surv(tt, follow) ~ 1, data = sample_tbl, conf.type = 'log-log')

ggsurvplot(follow_km)

follow_km %>% summary()

## Call: survfit(formula = Surv(tt, follow) ~ 1, data = sample_tbl, conf.type = "log-log")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5      4       1     0.75   0.217       0.1279        0.961
##     6      3       1     0.50   0.250       0.0578        0.845
##     7      1       1     0.00     NaN           NA           NA

2.2.1 The Hazard Estimator

To obtain more smooth functions we use kernel density estimator techniques.

sample_muhaz <- muhaz(sample_tbl$tt, sample_tbl$cens, max.time = 7)

plot(sample_muhaz)

broom has tidying methods for muhaz() and this allows us to create plots with ggplot2

muhaz_tidy_tbl <- sample_muhaz %>% tidy()

muhaz_tidy_tbl %>% glimpse()

## Observations: 101
## Variables: 2
## $ time     <dbl> 0.00, 0.07, 0.14, 0.21, 0.28, 0.35, 0.42, 0.49, 0.56, 0.63...
## $ estimate <dbl> 0.000000000, 0.000000000, 0.000000000, 0.000871904, 0.0034...

We have estimates of the hazard function now and so can plot it.

ggplot(muhaz_tidy_tbl) +
    geom_line(aes(x = time, y = estimate)) +
    expand_limits(x = 0, y = 0) +
    xlab("Time") +
    ylab("Estimated Hazard")

2.2.2 Boundary Corrections

By default, muhaz() corrects for the boundary on both sides, but we may not wish this. To get estimates without this correction, we add it as an argument to the function call.

sample_nocorr_muhaz <- muhaz(sample_tbl$tt, sample_tbl$cens, max.time = 7
                            ,b.cor = 'none')

plot(sample_nocorr_muhaz)

sample_nocorr_muhaz_tidy_tbl <- sample_nocorr_muhaz %>% tidy()

ggplot(sample_nocorr_muhaz_tidy_tbl) +
    geom_line(aes(x = time, y = estimate)) +
    expand_limits(x = 0, y = 0) +
    xlab("Time") +
    ylab("Estimated Hazard")

To help ensure these smoothed estimates are capturing the correct aspects of the data, we also have the equivalent pehaz() functions - giving histogram estimates of the hazards, much how histogram and kernel estimates are discrete and continuous analogies of one another.

2.3.1 No Boundary Correction

We also want to compare the results without the boundary correction.

muhaz_nocorr_surv_tbl <- sample_nocorr_muhaz_tidy_tbl %>%
    mutate(dt  = c(0, diff(time))
          ,S_t = exp(-cumsum(estimate * dt))
    )

ggsurvplot(example_km)$plot +
    geom_line(aes(x = time, y = S_t), colour = 'blue', data = muhaz_surv_tbl) +
    geom_line(aes(x = time, y = S_t), data = muhaz_nocorr_surv_tbl)

We see that removing the boundary corrections can cause big discrepancies between the smoothed and discrete estimates (at least in this case)

2.4.1 Plotting the Curves

Having created the KM estimates, we first plot the curves to check for differences.

telco_all_km   %>% ggsurvplot(censor = FALSE)

telco_vmail_km %>% ggsurvplot(censor = FALSE)

telco_intl_km  %>% ggsurvplot(censor = FALSE)

The estimators comes with upper and lower confidence bounds, and it may be useful to include these estimates when determining the differences. We will look at a more formal way of discriminating between curves in a while.

telco_all_km   %>% ggsurvplot(censor = FALSE, conf.int = TRUE)

telco_vmail_km %>% ggsurvplot(censor = FALSE, conf.int = TRUE)

telco_intl_km  %>% ggsurvplot(censor = FALSE, conf.int = TRUE)

2.4.2 Combining KM Estimators

This is useful, but it would be better to combine all these together. To do this, we need the tables of the estimator data, and stack them together.

Fortunately, the broom package is compatible with basic survival analysis routines, so this provides us with some utilities to help with this.

If we need to access the KM data directly, we can use tidy() to get a look:

telco_all_km_data   <- telco_all_km   %>% tidy()
telco_all_km_data   %>% glimpse()

## Observations: 218
## Variables: 8
## $ time      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
## $ n.risk    <dbl> 5000, 4989, 4987, 4979, 4976, 4974, 4972, 4967, 4965, 496...
## $ n.event   <dbl> 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 2, 0, 1, ...
## $ n.censor  <dbl> 10, 1, 7, 2, 2, 2, 5, 2, 3, 3, 6, 4, 9, 1, 5, 9, 6, 4, 8,...
## $ estimate  <dbl> 0.999800, 0.999600, 0.999399, 0.999198, 0.999198, 0.99919...
## $ std.error <dbl> 0.000200020, 0.000283183, 0.000347001, 0.000400944, 0.000...
## $ conf.high <dbl> 1.000000, 1.000000, 1.000000, 0.999984, 0.999984, 0.99998...
## $ conf.low  <dbl> 0.999408, 0.999045, 0.998720, 0.998414, 0.998414, 0.99841...

telco_vmail_km_data <- telco_vmail_km %>% tidy()
telco_vmail_km_data %>% glimpse()

## Observations: 404
## Variables: 9
## $ time      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 14, 15, 16, 17, 1...
## $ n.risk    <dbl> 3677, 3670, 3669, 3664, 3662, 3660, 3658, 3655, 3654, 365...
## $ n.event   <dbl> 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, ...
## $ n.censor  <dbl> 6, 0, 4, 1, 2, 2, 3, 1, 3, 2, 4, 6, 1, 5, 9, 4, 4, 5, 3, ...
## $ estimate  <dbl> 0.999728, 0.999456, 0.999183, 0.998911, 0.998911, 0.99891...
## $ std.error <dbl> 0.000271998, 0.000385030, 0.000471756, 0.000545035, 0.000...
## $ conf.high <dbl> 1.000000, 1.000000, 1.000000, 0.999978, 0.999978, 0.99997...
## $ conf.low  <dbl> 0.999195, 0.998702, 0.998260, 0.997844, 0.997844, 0.99784...
## $ strata    <chr> "vmailplan=no", "vmailplan=no", "vmailplan=no", "vmailpla...

telco_intl_km_data  <- telco_intl_km  %>% tidy()
telco_intl_km_data  %>% glimpse()

## Observations: 375
## Variables: 9
## $ time      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
## $ n.risk    <dbl> 4527, 4516, 4515, 4508, 4507, 4505, 4503, 4498, 4496, 449...
## $ n.event   <dbl> 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, ...
## $ n.censor  <dbl> 10, 1, 6, 1, 2, 2, 5, 2, 3, 3, 6, 3, 9, 1, 3, 8, 6, 4, 7,...
## $ estimate  <dbl> 0.999779, 0.999779, 0.999558, 0.999558, 0.999558, 0.99955...
## $ std.error <dbl> 0.000220921, 0.000220921, 0.000312845, 0.000312845, 0.000...
## $ conf.high <dbl> 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.00000...
## $ conf.low  <dbl> 0.999346, 0.999346, 0.998945, 0.998945, 0.998945, 0.99894...
## $ strata    <chr> "intlplan=no", "intlplan=no", "intlplan=no", "intlplan=no...

We stack the data and create a line plot of each of the survival curves.

stacked_km_tbl <- list(telco_all_km_data %>% mutate(strata = 'ALL')
                      ,telco_vmail_km_data
                      ,telco_intl_km_data
                       ) %>%
    bind_rows()


ggplot(stacked_km_tbl) +
    geom_line(aes(x = time, y = estimate,  colour = strata)) +
    geom_line(aes(x = time, y = conf.high, colour = strata), linetype = 'dashed') +
    geom_line(aes(x = time, y = conf.low,  colour = strata), linetype = 'dashed')

This chart is not as useful as it could be - it is too ‘busy’. We need to simplify.

Instead of including everything, we will look at the international plan curves along with the estimator calculated from the full dataset.

ggplot(stacked_km_tbl %>% filter(grepl('ALL|intl', strata))) +
    geom_line(aes(x = time, y = estimate,  colour = strata)) +
    geom_line(aes(x = time, y = conf.high, colour = strata), linetype = 'dashed') +
    geom_line(aes(x = time, y = conf.low,  colour = strata), linetype = 'dashed')

From this, it looks like people without international plans are similar to the population as a whole, but those with plans have a much worse churn effect. This type of information will be useful later.

We can do a similar thing for the voice-mail plans.

ggplot(stacked_km_tbl %>% filter(grepl('ALL|vmail', strata))) +
    geom_line(aes(x = time, y = estimate,  colour = strata)) +
    geom_line(aes(x = time, y = conf.high, colour = strata), linetype = 'dashed') +
    geom_line(aes(x = time, y = conf.low,  colour = strata), linetype = 'dashed')

This split is more mixed, but it appears that those with a voicemail plan have a better churn experience than the population, but those without are similar to the population as a whole.

3.1.1 Voicemail Differences

We start by comparing the voicemail plans.

telco_vmail_survdiff <- survdiff(Surv(accountdur, churned) ~ vmailplan
                                ,data = telco_churn_tbl)

telco_vmail_survdiff %>% print()

## Call:
## survdiff(formula = Surv(accountdur, churned) ~ vmailplan, data = telco_churn_tbl)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## vmailplan=no  3677      605      523      13.0      50.1
## vmailplan=yes 1323      102      184      36.8      50.1
## 
##  Chisq= 50.1  on 1 degrees of freedom, p= 1e-12

telco_vmail_survdiff %>% glance()

## # A tibble: 1 x 3
##   statistic    df  p.value
##       <dbl> <dbl>    <dbl>
## 1      50.1     1 1.48e-12

Such a small \(p\)-value suggests the differences in the curves are unlikely to be due to chance, and this gels with what we observed on the plots of the curves.

3.1.2 International Calls Differences

Now we compare the differences between plans with and without international calls.

telco_intl_survdiff <- survdiff(Surv(accountdur, churned) ~ intlplan
                               ,data = telco_churn_tbl)

telco_intl_survdiff %>% print()

## Call:
## survdiff(formula = Surv(accountdur, churned) ~ intlplan, data = telco_churn_tbl)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## intlplan=no  4527      508    638.6      26.7       277
## intlplan=yes  473      199     68.4     249.3       277
## 
##  Chisq= 277  on 1 degrees of freedom, p= <2e-16

telco_intl_survdiff %>% glance()

## # A tibble: 1 x 3
##   statistic    df p.value
##       <dbl> <dbl>   <dbl>
## 1      277.     1       0

3.3.1 Discretizing custservicecalls

Having decided to bin our calls, we now need to decide how we go about this. One way is to split the data by quantiles, so we can try this first.

Should this not prove satisfactory, it may help inform more suitable ways to bin this value.

telco_churn_tbl %>% count(custservicecalls) %>% mutate(cuml = cumsum(n))

## # A tibble: 10 x 3
##    custservicecalls     n  cuml
##               <dbl> <int> <int>
##  1                0  1023  1023
##  2                1  1786  2809
##  3                2  1127  3936
##  4                3   665  4601
##  5                4   252  4853
##  6                5    96  4949
##  7                6    34  4983
##  8                7    13  4996
##  9                8     2  4998
## 10                9     2  5000

We see this is a very skewed distribution to the lower end of call counts, so ntile() is not as effective here as it will not properly separate the counts. Instead we use cut() as that gives us more control over the binning.

telco_churn_cat_tbl <- telco_churn_tbl %>%
    mutate(calls_cat = cut(custservicecalls
                          ,breaks = c(0, 1, 2, 5, Inf)
                          ,right = FALSE)
    )

telco_callscat_km <- survfit(Surv(accountdur, churned) ~ calls_cat
                            ,data = telco_churn_cat_tbl)

telco_callscat_km %>% ggsurvplot(censor = FALSE, conf.int = TRUE)

It looks like there is not much difference for lower counts, but there is a worse churn experience when the number of customer service calls is 5 or higher.

telco_callscat_survdiff <- survdiff(Surv(accountdur, churned) ~ calls_cat
                                   ,data = telco_churn_cat_tbl)

telco_callscat_survdiff %>% print()

## Call:
## survdiff(formula = Surv(accountdur, churned) ~ calls_cat, data = telco_churn_cat_tbl)
## 
##                      N Observed Expected (O-E)^2/E (O-E)^2/V
## calls_cat=[0,1)   1023      121    145.1      3.99      5.04
## calls_cat=[1,2)   1786      190    254.7     16.45     25.82
## calls_cat=[2,5)   2044      306    285.4      1.49      2.51
## calls_cat=[5,Inf)  147       90     21.8    213.33    220.88
## 
##  Chisq= 236  on 3 degrees of freedom, p= <2e-16

telco_callscat_survdiff %>% glance()

## # A tibble: 1 x 3
##   statistic    df p.value
##       <dbl> <dbl>   <dbl>
## 1      236.     3       0

4.3.1 R-Squared

The standard measure of model fit is \(R^2\), the proportion of variance explained by the model. This does not quite work for survival models, but we can use a modified \(R^2\) statistic, defined for Cox-PH models as follows:

\[ R^2 = 1 - \left( \frac{l(0)}{l(\hat{\beta})} \right)^{\frac{2}{n}} \]

where \(l(\beta)\) is the log-likelihood function.

4.3.2 Concordance Index

The concordance index, or \(c\)-index, or \(c\)-statistic is the proportion of pairwise comparisons correctly identified by the model. To ensure a fair comparison at least one of the rows must not be censored.

We compare comparison times and check that our model correctly predicted the order of times between the compared datapoints - a correct value adds 1 to the total, an incorrect one adds 0. We then divide the total by the number of comparisons made to get our statistic.

Thus, a concordance index of 0.5 is equivalent to a random guess and so is our baseline value. A good survival model will have a \(c\)-index between 0.6 and 0.7 - anything higher is unlikely due to the censored nature of our data.

4.3.3 Assessing Our First Model

So how good is our model?

telco_model01_coxph %>% summary()

## Call:
## coxph(formula = Surv(accountdur, churned) ~ vmailplan, data = telco_churn_tbl)
## 
##   n= 5000, number of events= 707 
## 
##                coef exp(coef) se(coef)     z Pr(>|z|)    
## vmailplanyes -0.741     0.477    0.107 -6.92  4.7e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## vmailplanyes     0.477        2.1     0.387     0.588
## 
## Concordance= 0.565  (se = 0.008 )
## Rsquare= 0.011   (max possible= 0.88 )
## Likelihood ratio test= 56.6  on 1 df,   p=5e-14
## Wald test            = 47.8  on 1 df,   p=5e-12
## Score (logrank) test = 50  on 1 df,   p=2e-12

telco_model01_coxph %>%
    glance() %>%
    select(n, nevent, r.squared, concordance, std.error.concordance)

## # A tibble: 1 x 5
##       n nevent r.squared concordance std.error.concordance
##   <int>  <dbl>     <dbl>       <dbl>                 <dbl>
## 1  5000    707    0.0113       0.565               0.00757

telco_model01_coxph %>% tidy()

## # A tibble: 1 x 7
##   term         estimate std.error statistic  p.value conf.low conf.high
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
## 1 vmailplanyes   -0.741     0.107     -6.92 4.67e-12   -0.951    -0.531

This simple model with one binary variable gives us a \(c\)-index of 0.565. The \(R^2\) is essentially zero (0.011), and we will check how both values correlate as we iterate upon the model-building.

4.5.1 Fitting the Null Model

We first work from a null CoxPH model that we fit - the ‘null’ model:

telco_model00_coxph <- coxph(Surv(accountdur, churned) ~ 1
                            ,data = telco_churn_tbl
                             )

telco_model00_coxph %>% summary()

## Call:  coxph(formula = Surv(accountdur, churned) ~ 1, data = telco_churn_tbl)
## 
## Null model
##   log likelihood= -5293.32 
##   n= 5000

4.5.2 Martingale Residuals

We define our basic residual \(m_i\) as follows:

\[ m_i = \delta_i - \hat{H}_0(t_i) \exp(x_i \hat{\beta}) \]

The residual is the difference between the observed value of the whether the event occurred and its expected value under the Cox model.

telco_model00_resid_m <- telco_model00_coxph %>% residuals(type = 'martingale')

telco_model00_resid_m %>% summary()

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.2069 -0.1805 -0.0704  0.0000 -0.0166  0.9998

We have zero-mean, but there is skew in the data as the values are bounded above at 1, but unbounded below. The easiest way to check the shape is look at a histogram of the values.

ggplot() +
    geom_histogram(aes(x = telco_model00_resid_m), bins = 50) +
    xlab('Residual Value') +
    ylab('Count') +
    ggtitle('Histogram of Null-Model Martingale Residuals')

ggplot() +
    geom_histogram(aes(x = telco_model00_resid_m^2), bins = 50) +
    xlab('Squared Residual Value') +
    ylab('Count') +
    ggtitle('Histogram of Null-Model Squared Martingale Residuals')

4.5.3 Deviance Residuals

Martingale residuals cannot be used in sums of squares methods to assess goodness-of-fit - instead we define another quantity when the model is correct is symmetrically distributed with an expectation of zero. These ‘deviance’ residuals, \(d_i\) are defined as follows:

\[ d_i = \text{sign}(m_i) {-2 [m_i + \delta_i \log (\delta_i - m_i)]}^{1/2} \]

telco_model00_resid_d <- telco_model00_coxph %>% residuals(type = 'deviance')

telco_model00_resid_d %>% summary()

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -1.554  -0.601  -0.375  -0.159  -0.182   3.877

The mean is small, but non-zero. We inspect the distribution of values to see what it looks like.

ggplot() +
    geom_histogram(aes(x = telco_model00_resid_d), bins = 50) +
    xlab('Residual Value') +
    ylab('Count') +
    ggtitle('Histogram of Null-Model Deviance Residuals')

ggplot() +
    geom_histogram(aes(x = telco_model00_resid_d^2), bins = 50) +
    xlab('Squared Residual Value') +
    ylab('Count') +
    ggtitle('Histogram of Null-Model Squared Deviance Residuals')

Neither look distributed as we would like, suggesting this model is not good, but we will return to this topic later.

4.5.4 Residual Plots

Much like for linear models, these residuals are useful to help decide which variables may have predictive power in our survival model. We do this by plotting the values of the residual against the variable values.

To do this, we first add the residuals to the original dataset so that we can construct some of the residual plots.

telco_null_assess_tbl <- telco_churn_tbl %>%
    mutate(resid_m = telco_model00_resid_m
          ,resid_d = telco_model00_resid_d
           )

Now that we have put the table together we can construct some two-way plots of the variable against the martingale residual:

ggplot(telco_null_assess_tbl) +
    geom_boxplot(aes(x = vmailplan, y = resid_m)) +
    xlab('vmailplan') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against vmailplan')

ggplot(telco_null_assess_tbl) +
    geom_boxplot(aes(x = intlplan, y = resid_m)) +
    xlab('intlplan') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against intlplan')

ggplot(telco_null_assess_tbl) +
    geom_boxplot(aes(x = areacode, y = resid_m)) +
    xlab('areacode') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against areacode')

ggplot(telco_null_assess_tbl) +
    geom_boxplot(aes(x = custservicecalls, y = resid_m, group = custservicecalls)) +
    xlab('custservicecalls') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against custservicecalls')

custservicecalls is a numeric variable, but it does not take many values so we can create boxplots rather than a more typical scatterplot.

For continuous variables, we look at those scatterplots:

ggplot(telco_null_assess_tbl) +
    geom_point(aes(x = vmailmessage, y = resid_m), size = 1, alpha = 0.1) +
    xlab('vmailmessage') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against vmailmessage')

ggplot(telco_null_assess_tbl) +
    geom_point(aes(x = daymins, y = resid_m), size = 1, alpha = 0.1) +
    xlab('daymins') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against daymins')

ggplot(telco_null_assess_tbl) +
    geom_point(aes(x = nightmins, y = resid_m), size = 1, alpha = 0.1) +
    xlab('nightmins') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against nightmins')

ggplot(telco_null_assess_tbl) +
    geom_point(aes(x = intlmins, y = resid_m), size = 1, alpha = 0.1) +
    xlab('intlmins') +
    ylab('Martingale Residual') +
    ggtitle('Residual Plot of Martingal Residual against intlmins')

No strong patterns stand out here, but that does not necessarily rule them out of candidate models. We will not know the predictive power till we try them in a model.

4.6.1 Residual Plots

Tests for the proportional hazards assumption reduce to checking if the parameters from the model are constant over time. This task is made easier by the construction of ‘Schoenfeld’ residuals, which we will not discuss here, but are easily calculated in R by using the cox.zph() function.

telco_model01_coxzph <- telco_model01_coxph %>% cox.zph(transform = 'km')

telco_model01_coxzph %>% print()

##                 rho chisq     p
## vmailplanyes 0.0597  2.52 0.112

telco_model01_coxzph %>% ggcoxzph()

The plot and \(p\)-value suggest the assumptions are not horrible here.

We now move on to the second model we built.

telco_model02_coxzph <- telco_model02_coxph %>% cox.zph(transform = 'km')

telco_model02_coxzph %>% print()

##                rho chisq     p
## intlplanyes 0.0269 0.512 0.474

telco_model02_coxzph %>% ggcoxzph()

The assumption also seems sound here.

We now try the combined model, model03:

telco_model03_coxzph <- telco_model03_coxph %>% cox.zph(transform = 'km')

telco_model03_coxzph %>% print()

##                 rho chisq     p
## vmailplanyes 0.0547 2.100 0.147
## intlplanyes  0.0298 0.626 0.429
## GLOBAL           NA 2.780 0.249

telco_model03_coxzph %>% ggcoxzph()